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ABSTRACT 

We present a novel solver for an analogue to Poisson's equation in the framework of 
modified Newtonian dynamics (MOND). This equation is highly non-linear and hence 
standard codes based upon tree structures and/or FFT's in general are not applicable; 
one needs to defer to multi-grid relaxation techniques. After a detailed description of 
the necessary modifications to the cosmological iV-body code AMIGA (formerly known 
as MLAPM) we utilize the new code to revisit the issue of cosmic structure formation 
under MOND. We find that the proper (numerical) integration of a MONDian Pois- 
son's equation has some noticable effects on the final results when compared against 
simulations of the same kind but based upon rather ad-hoc assumptions about the 
properties of the MONDian force field. Namely, we find that the large-scale structure 
evolution is faster in our revised MOND model leading to an even stronger clustering 
of galaxies, especially when compared to the standard ACDM paradigm. 

Key words: galaxy: formation - methods: A^-body simulations - cosmology: theory 
- dark matter - large scale structure of Universe 



1 INTRODUCTION 

Modified Newt onian dynamics (MOND) was proposed by 
iMilgronj l|l983l ) as an alternative to Newtonian gravity to 
explain galactic dynamics without the need for dark matter. 
Althougli current cosmological observations point to the ex- 
istence of vast amounts of non-bary onic dark matter in the 
Universe (e.g. iKomatsu et ahlboOSl ). it remains interesting 
to explore other alternatives, especially as not all of the fea- 
tures of CDM models appear to match observational data 
(e.g., the "missin g satellite problem" (|Klvpin et al.l Il999l : 
Moore et al.lll999f) a nd the so-called "cu sp-core crisis" (e.j 



de Blok et al.l l2003l : ISwaters et al] 120031 )). In that regards 



it appears important to look for tests able to discriminate 
between MOND and Newtonian gravity, especially in the 
context of cosmology now that there exist various relativis- 
ticformulation of the MOND theory l|Bekensteinll2004 IZhad 
I2OO7I . I2OO8I ). However, progress in that field has been ham- 
pered by the fact that MOND is a non-linear theory making 
any analytical predictions as well as numerical simulations 
a tedious task. To date, only a few codes exist that actually 
solve the MONDian a nalogue to Poisson equation in a non- 
cosmolo gical conte xt (Brad a fc MilgromI 1 19991 : iNipoti et al.l 
[2007; Ti ret fc Com bes 200 1), none of which is publically 
available. We are augmenting this list by making available 
a new solver for the MONDian Poisson's equation primar- 
ily designed to work in a cosmological context but readily 
adjusted to allow for simulations of isolated galaxies. 



Until recently MOND was merely a heuristic theory tai- 
lored to fit rotation curves with little (if any) predictive 
power for cosmological structure formation. One of the most 
severe problems for the general appreciation and acknowl- 
edgment of MOND as a "real" theory (and a conceivable 
replacement for dark matter) was the lack of success to for- 
mulate the theory in a general relativistic manner. This sit- 
uation though changed during the last couple of years and 
at present there are a number of covariant theories (e.g. 
(|Bekensteinll2004 IZhaoll2007l . [ 20081 )) whose non-relativistic 
weak acceleration limit accords with MOND while its non- 
relativistic strong acceleration regime is Newtonian. How- 
ever, there remains a lot to be done with regards to the rel- 
ativistic formulations of MOND in order to a get completly 
acceptable theory. But nevertheless has MOND reached a 
stage of development in wich we can think in doing cos- 
mology with reducing the need of unjustifiable assumptions. 
Furthermore, the study of cosmological structure formation 
in the non-linear regime will provide new constrains on such 
generalizations of MOND. 

The first valiant attempts at simulating cosmic struc- 
ture formatio n un der the influence of MOND were done by 
iNusseij (|2002l ) and lKnebe et all l|2004 KG04 from now on). 
While their studies provided great insights into the non- 
linear clustering behaviour of MONDian objects (we dare to 
call them galaxies as none of these simulations included the 
physics of baryons) they were still based upon some rather 
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unjustifyable assumptions. The main objective of this paper 
now is to refine the implementation of MOND in theiV- 
body code AMIGA (successor of MLAPM. iKnebe et alj (j200ll l'l. 
i.e. we modified the code to numerically integrate a MON- 
Dian analogue to Poisson's equation. This equation is a 
highly non-linear partial differential equation whose solu- 
tion is non-trivial to obtain. Only sophisticated multi-grid 
relaxation techniques (on adaptive meshes of arbitrary ge- 
ometry) are capable of tackling this task. We argue that only 
when MOND has been thoroughly and "properly" studied 
and tested against ACDM can we safely either rule it out for 
once and always or confirm this rather venturesome theory. 
The theory has become a valid competitor to dark mat- 
ter and it therefore only appears natural - if not manda- 
tory - to (re-)consider its implications. We developed a tool 
and the subsequently necessary analysis apparatus allow- 
ing to test and discriminate cosmological structure forma- 
tion in a MONDian Universe from the standard dark matter 
paradigm and used it to study a particular model of MOND 
theory. 



2 THE MONDIAN EQUATIONS 

Despite the recent progress made in obtaining, refining and 
studying covariant formulati o ns of the MOND theory (cf. 
iBekensteinlliooi : IZhaol [2OO7I . [ 20081 ) . there are still certain 
ambiguities in contriving a closed a set of equations-of- 
motion suitable for an A^-body code. We will work out such 
a set in this Section. 

Non-MONDian Cosmology We like to remind the 
reader that in a cosmological A^'-body code one integrates 
the comoving equations of motion 

ma'' 

p — — mV$jv 

which are completed by Poisson's equation for the Newto- 
nian potential responsible for peculiar accelerations 

V-[V$iv(x)] = l^(p(x)-p) . (2) 
a 

In these equations a is the cosmic expansion factor, x = r/a 
the comoving coordinate, p the canonical momentum, V- the 
divergence operator with respect to x, p (p) the comoving 
(background) density, and 



givW = -V-I'iv(x) 



(3) 



the Newtonian peculiar acceleration field in comoving coor- 
dinates. 

The next step now is to define an analogon to Eq. ((5)) 
that takes into account the effects of MOND in cosmology. 

MONDian Cosmology In order to get an equation for the 
comoving peculiar MONDian potential $m we have to make 
a decision about which covariant generalization of MOND 
we want to use. Assuming the p function to be constant 
in space, a non-relativistic limit of the covariant theory by 
IZhad (|2008t ) can be written as follows 



where fi{x) is p ^ 1 for a; 2> 1 (Newt onian limit) an d 
p ^ X for a; < 1 (MONDian limit) (cf. iMilgromI Il983l ) . 
We further took the liberty to encode the MONDian ac- 
celeration scale 7(a) as a (possible) function of the cos- 
mic expansion factor a for reasons that will become clear 
later on (cf. Section Q. The most naive choice would be 
7(a) ~ go = 1.2 X 10~*cm/sec^ whereas other theori e s ma; 
lead to different dependencies; for instance, in IZhaol (|200 ' 
7(a) is given as 7(a) = a^^^go. 

In analogy to Eq. |[3]) we define 



gM(x) 



-V'l>Al(x) 



(5) 



as the MONDian peculiar acceleration field in comoving co- 
ordinates. 

The Curl-Field The relation between the MONDian and 
t he Newtonian force, i.e. Eq . ([HJ and Eq. ([S}, is given by 
(|Bekenstein fc Milgro"mlll984 ) 



go 



Sm =Sn + C 



where an otherwise unspecified curl-field 
C = V X h , 



(6) 



(7) 



1V£m| 
aj{a) 



V'I'j 



4ttG 



ip-p) 



(4) 



appears that has been shown to vani sh for any kind of sym- 
metry (jBekenstein fc MilgromI Il984l ). Previous studies on 
cosmological structure formation under MOND neglected 
the curl-field C as this allows to use a standard solver for the 
Newtonian Poisson's equation ([2} and then an inversion of 
Eq. (O to obtain the MONDian forces (.Nusser. .200a. KG04). 
However, this leads to parallel Newtonian and MONDian 
forces which is not necessarily the case! The more correct 
approach is to directly solve Eq. Q and in the following 
Section we present a novel solver that numerically solves for 
the MONDian potential "1>m to be used with Eq. ((5} in or- 
der to integrate the equations-of-motion ([l| for N particles. 
The solution to Eq. Q can in turn be used together with 
Eq. ([6} to actually study the curl-field C and its effects. 



3 THE CODE 

Conventional Poisson solvers used throughout (computa- 
tional) cosmology are no longer applicable to Eq. Q. All 
of the standard methods such as Fast-Fourier- Transform 
(FFT) based Particle Mesh (PM) codes, tree codes, ex- 
pansion codes, and their variants rely on the linearity of 
Poisson's equation. The corresponding MONDian Poisson 
equation (|4]) though is a non-linear partial differential equa- 
tion which requires more subtle and refined techniques to be 
solved numerically. We are now going to elaborate upon the 
steps required to adjust a given relaxation solver to account 
for the complexity of the MONDian Poisson's equation and 
any arbitray MONDian interpolation function p{x). 

3.1 Solving the MONDian Poisson's equation via 
multi-grid relaxation 

The code used is a modification of the open sou rce cosmolog- 
ical c ode AMIGA that is the sucessor of MLAPM (|Knebe et al.l 
l200lf ). We conserve the equations of motion of the origi- 
nal code (i.e. Eq. ([T}) but we now numerically integrate the 
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MONDian Poisson's equation Q. The fact that the new 
equation for the potential is non-linear prevents us from us- 
ing the standard methods (cf. above) ; we rather need to ad- 
here to a multi-grid relaxation techn ique (e.g.. iBrandt 197?! : 
iPress et al.lll99^ : iKnebe et al.ll20oil ). The method consists 
of discretizing the equation on a grid and solving the non 
linear system of algebraic equations by relaxing a trial po- 
tential "I> un til convergence. We adhere to the discretisation 
proposed bv iBrada fc MilgromI l|l999t ) and already used by 
iTiret fc Combed l|2007l ) 



(/Xi+l/2(<^i+l - 4>) + Mi-l/2(0 - (Pi-l) + 
Hj+l/2{<l>j+\ -<!})+ Mj-l/2(<^ - </>j-l) + 

^fc+i/2(<^fc+i - 4>) + ^^k-\/2{<i> - <t}k-i))/h? 

4ttG 



(8) 



where Pij^k is the density contrast in cell {i,j,k) on a grid 
with spacing h, ij> is the potential in cell {i,j,k) with <l)i±i 
the potential in the neighbouring cells (i± l,j, k). The same 
terminology applies to the other two spatial dimensions. The 
coefficient Hi+i/i is the evaluation of the ^i(x) function at 
the point {i + 1/2, j, k), etc. 
Reordering Eq. ((Sjl we get 

pi+\/24>i+l + Mi-l/20i-l + 
Mj+l/2<^j + l + Mj-l/20j-l + 
Mfc+l/2<^fc + l + P'k-l/2<t>k-l — 
iPi+1/2 + Mi-1/2 + Mj + 1/2 + 
Mj-1/2 + Mfc+1/2 + Pk^l/2)4> = 

47rG 



-Pi,j,k 



(9) 



where we can see that this way of discretising the equation 
is similar to the standard discretization of the Newtonian 
Poisson's equation, but the potential is "weighted" by the /i 
function. Given the nature of the ^ function under MOND 
it can be set to unity for accelerations larger than go and 
we then recover the standard discretisatio n of the Newtonia n 
Poisson's equation (cf. equation (11) in Knebe et al.l l|200lh ). 

The gradient of the potential, i.e. Vi^, in the argument 
of i-i{x) needs to be disctretized, too. We choose the following 
form 



((V<^).),+l/2 - ^ 



(10) 



/2 



+ 1/2 



(<^i+l,j+l + (jti + l) — (0^+l,j-l + <t>j~l) 
Ih ' 

(<^i+i,fc+i + 0fc+i) — {(l>i+i,k-i + 4>k-i) 

The nodes used in this discretisation are highlihted in Fig- 
ure [T] 

There are two possible choices for the relaxation proce- 
dure. We can either re-order terms in Eq. ((9)1 again (while 
freezing the coefficients fJ.i±i/2, Pj±i/2, Pk±i/2) in a way that 
we get 

<t) = /(</>i±l,0j±l,<?!'fc±l,Mi±l/2,Mj±l/2,A*fc±l/2,Pij,fc), (11) 

or we can apply one step of a Newton- 
Raphson root-finding algorithm to — ~ 




Figure 1. Stencil (27 points) used in the discretization of the 
MONDian Poisson's equation (cf. Eq. Jojl). The argument of the 
/j(x) function is approximated in points indicated by the open 
circle, i.e. in-between the standard nodes. The triangles show the 
cells used to calculate the gradient in that point (cf. Eq. ((TOj). 



f ((l>i±l, 07±l i 0fc±li A'»+l/2i Mi + l/2i ^■k + l/2, Pi.i.k) 

([Press et al.l [X992). These two alternatives show simi- 
lar rate of convergence that goes in favour of one or the 
other according to the density. Our method of choice is the 
iterative procedure as it involves fewer calculations. 

We still need to define a "colouring" -scheme, i.e. a way 
of how to sweep through the whole grid and (iteratively) 
update the potential <^ in a given cell k). The ordering 
of the iterations is given by a generalization of the standard 
two-colour/red- black method and sketched in Figure [2] One 
iteration step is complete after sweeping eight times trought 
the grid updating those nodes per sweep that have identical 
numbers as indicated in Figure (2] 

We further coded additional iteration schemes (e.g., 
lexicographic and zebra in three directions) and the code 
switches automatically between them in extreme cases for 
which the convergence with the 8 colours is slow. For 
an elaborate discussion of multi-grid relaxation techniques 
and colo uring-schemes in particular we refer the interested 
reader to IWesselin"3 (|l992l ) . 

A numerical solution has converged after n itera- 
tions on a grid k if the norm 1 1 • 1 1 (mean or maximun value 
in the grid) of the residual 



AtvG 



ip-p) 



(12) 



is small compared against the norm of the truncation error 

r*' = L'=-i(7?$J;)-i?L''($^) (13) 

i.e. Ile'^ll < 0.1||r'''||. The value 0.1 has been determined 
heuristically. Here, denotes the discretization of the left 
hand side of Eq. (|4]) on a grid k, L*'"^ the same discretization 
in the next coarser grid and R is the restriction operator used 
to interpolate values from the grid k to the grid k ~ 1. 
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Figure 2. Eight-colour scheme used for obtaining an iterative 
solution of Eq. Q (cf. also Eq. mil l'). The numbers indicate the 
ordering by which the nodes that are updated. 



3.2 Tests 

In order to verify and gauge the credibility of our numerical 
solver we performed two complementary tests. The first sim- 
ply assesses the recovery of the potential for a given solution 
of the MONDian Poisson's equation. The second test checks 
whether or not we recover the same temporal evolution of a 
(cosmological) simulation when setting the MONDian accel- 
eration go to such a low value that it will not affect structure 
formation (Newtonian limit). 

3.2.1 Static test 

As a first (static) test for the numerical potential solver, we 
present the recovery of the analytical potential for a known 
analytical solution of the MONDian Poisson's equation @. 
To this extent we actually start with an analytical descrip- 
tion for a sp herically symmet ric potential akin to the Hern- 
quist profile (|Hernauistlll990l) 



GM 



- + y/GMgo In ( 



ro 



(14) 



We added a logarithm to the potential in order to mach the 
properties of the solution of the MOND equation (far from 
any source, the MONDian forces are proportional to In (r)). 
We then derive the density by substituting Eq. (|14p into 
Eq. (Hi 



P = -3 



+ 



go C go C • gl 



A 



B = 



GM 



VGM^ 

(r + rc)^ (r + rc) 

2GM ^/GMgo 

(r -I- Tc)-* (r + r-c)^ 

A 



(15) 



C = 1 + — 

90 

The test is performed with non-periodic boundary condi- 
tions on a regular 64"^ grid by fixing the solution on the 
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Figure 3. Relative error in the potential (upper panel) and the 
forces (lower panel) recovered with the solver on a 64'^ grid for 
a Hernquist profile as a function of radius. The solid (red) line 
is the diagonal across the box while the dots represent a random 
sample of all grid points. 



border to the analytical values. Figure|3]now shows the frac- 
tional error in the potential (upper panel) and force (lower 
panel). The constants in the potential for this particular test 
are: M = 10^°Mq, = 2.5 x 10"^ Mpc and ro = 0.01 Mpc. 
The box used was 0.01 Mpc. With this choice for the param- 
eters, the force gradually changes from a purely Newtonian 
regime in the central parts to deep MONDian in the outer 
regions. The fractional errors are calculated as the difference 
between the analytical potential and forces on the grid and 
the numerical solution obtained by our solver. The solid red 
line shows the error as a function of radius along the diago- 
nal in the box while the dots represent a random sample of 
all grid points^ We can clearly see that the error is never 
larger than 1 per cent indicating the excellent quality of our 
numerical integration scheme. 



^ Note that running along the 3-dimensional diagonal will coin- 
cide with the largest possible error and hence the solid line in the 
upper panel of Figure |3] marks the upper boundary of the error 
in the potential. 
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Figure 4. Density distribution at redsliift ^ = for the 64^ par- lOO looo loooo 

ticle test simulation in the Newtonian limit (30 = lO"-"^^ cm/s^, 
upper panel) and the corresponding Newtonian run (lower panel) . 

Figure 6. Mass function of the 64"^ particle test simulation at 
redshift 2 = 0. 



3.2.2 Dynamic test 

Having established the credibility of our numerical integra- 
tor for Eq. Q, we now test the temporal evolution of our 
code by using a cosmological simulation in the Newtonian 
limit, i.e. lowering the MONDian acceleration go to such a 
low value that it will not impact upon structure formation. 
For that purpose we used go = 10~^^cm/s^ and compare the 
final output to a simulation run with the standard Newto- 
nian A'^-body code. The simulation fascilitated 64'^ particles 
in a box of comoving side length 15h~^ Mpc. The cosmology 
is characterized by r2o=0.3, f2A,o=0.7, as=0.9, and h = 0.7 
but actually of no relevance for this particular test. 

A first visual impression of the density field in given 
in Figure |4] where we show the grey-scaled density at each 
particle position of the simulation at redshift z = 0. The 
differences are at best marginal indicating that the modified 
solver performs correctly in the Newtonian limit. A more 



quantitative comparison can be found in Figure [5] where we 
show the matter power spectrum at redshift z — 0. The 
curves are indistinguishable which is proven by the ratio of 
the two curves plotted as a dotted line (lowered for clar- 
ity by the factor 0.001). We further ran the MPI enabled 
halo finder AHlQ (KnoUmann et al. 2008, in prep.) over both 
simulations. The resulting mass functions of identified ob- 
jects is given in Figure |6] Again, there are at best differ- 
ences that are readily ascribed to numerical errors during 
the time integration due to our vastly different schemes of 
solving Poisson's equation in the two test runs. 



^ AHF is the successor of MHF introduced in I Gill et al. 
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Table 1. Model parameters. In all cases a value for the Hubble 
parameter of h = 0.7 was employed. For model OCBMond wc 
explicitly assumed C = whereas OCBMond2 is based upon 
a numerical integration of the MONDian Poisson's equation as 
described in Section |3] 

label f^o \) cr|"" o-8°™ 90 [cm/s^] 



ACDM 

OCBMond 

OCBMond2 



0.30 
0.04 
0.04 



0.04 
0.04 
0.04 



0.7 
0.0 
0.0 



0.88 
0.92 
0.92 



0.88 
0.40 
0.40 



1.2 xlO"* 
1.2 xlO"** 



4 COSMOLOGICAL SIMULATIONS 

To date there is only a limited number of cosmological simu- 
lations aiming at s tudying the e ffect sof MOND in a full cos- 
mological context (|Nusseij2002l . KG04). However, both these 
investigations "simply" solved for the Newtonian forces and 
modified them according to a non-cosmological MOND pre- 
scription prior to being used with the equations-of-motion 
(cf. Eq. (HI). They both showed that, under their respec- 
tive heuristic formulations, MOND can lead to large-scale 
clustering patterns aking to a ACDM model. 

In the present study we are going to abandon (at least) 
one of the assumptions made by these authors, namely that 
C = and hence that Eq. ^ can be inverted to obtain a 
function gM(gjv). Here we numerically integrate the MON- 
Dian analogue to Poisson's equation Q which directly leads 
to considering the influence of the yet unkown and unre- 
garded curl-field C in the process of structure formation. 
The "only" assumption we adhere to is that MOND does 
not affects fluct uations and leaves the backgrou nd cosmol- 
ogy intact (e.g.. lNusse3l2002l : iKnebe et al.ll2004 ). We defer 
to a later study that will make use of modifled Friedmann 
equations accounting for the effects of MOND. 



4.1 The Simulation Details 

Following KG04 and using th e input power spectra de - 
rived with the CMBFAST code (|Seliak fc Zalda rriaga"l996") 
we displace 128'^ particles from their initial positions 
on a regular lat t ice u sing the Zel'dovich approximation 
lEfstathiou et al.l 119851 ). The box size was chosen to be 
32h~^ Mpc on a side. This choice guarantees proper treat- 
ment of the fundamental mode which will still be in the 
linear regime at z — (cf. the scale turning non-linear at 
z = is roughly 20h~^ Mpc for the models under inves- 
tigation). The particulars (and terminology) of the models 
under investigation are summarized in Table [T] 

As we used the (Newtonian) Zeldovich approximation 
to generate the initial conditions we need to make sure that 
the universe is still Newtonian at the starting redshift. To 
this extent our choice for the (free) function 7(a) in Eq. (|4]) 
is 7(a) ~ ago- This ensures that /i(x) — 1 and hence Eq. (|4]) 
reduces to the Newtonian case Eq. ([5}. We need to acknowl- 
edge that this differs from the treatment in KG04 who ac- 
tually used 7(a) — go; their simulations did not start early 
enough and were in parts already MONDian at the initial 
redshift. However, as we will see later on during the analysis 
of the runs that our results hardly differ from the findings 
of KG04. We further like to mention that the present study 
does not aim at distinguishing one MOND theory from an- 
other! We merely apply our novel solver to one particular 



of the runs is rrip = 1.30 x 10^ h 



MOND model of our choice showing that this implementa- 
tion can lead to results that are at least comparable to the 
favoured ACDM structure formation scenario. Besides, we 
also like to study the effects of the yet unregarded and un- 
kown curl-field C upon gravitational structure formation, 
irrespective of the applicability of our particular MOND 
model to reality; we leave a detailed study of various MOND 
realisations to a future paper. 

The particles were evolved from redshift z — 50 until 
z — and in all three runs a force resolution of 6h~^ kpc was 
reached in the highest density regions. The mass resolution 

1 Mo for the ACDM model 
'^h~^ Mq for the two low-fio models, 
respectively. We output 26 snapshots of the particle posi- 
tions and velocities equally spaced in time t inbetween red- 
shifts z = 5 and z = 0. These outputs are then analysed 
with respect to their large-scale clustering patterns as well 
as properties of individual objects employing the MPI en- 
abled halo finder AHF again. For that purpose we though 
switched off the unbinding procedure; we rather collect all 
particles within a certain region of a local density peak and 
refer to them as "object particles" defining the properties of 
that object. To be more precise, AHF locates density max- 
ima in the simulation by invoking the original adaptive-mesh 
hierarchy again used during the simulation procedure. For 
each of these peaks we step out in logarithmically spaced 
radial bins until the mean density inside that (spherical) 
sphere drops below a fiducial value A x p6 where pb is the 
cosmological background density, i.e. 



M(< R) 



= Ax Pb, 



(16) 



defines the edge R of an object. However, one needs to care- 
fully choose the correct overdensity A which is much higher 
for the OCBMond models due to the low Oq value. The 
parameters used a re A = 340 for ACDM and A = 2200 
for OCBMond (see ICros^l 19971 . Appendix C, and references 
therein). 



4.2 The Matter Field 

We start with inspecting the matter density field of all three 
runs in Figure [T] There we can see that at redshift z — 
the differences amongst the models are rather small. How- 
ever, there are major differences at larger redhifts due to the 
fact that the formation of structures under MOND kicks in 
later and then at an increased growth rate (i.e. 5 cx a'^ in 
MOND as op p osed to 5 oc a for Newtonian ph y sics, e.g. 
S andersI l|200ll ) ; iNussej hOOt ; iKnebe et al.1 (|2004| ) ; ISanderj 
( 120081 )'). But the more interesting observation is the compar- 
ison between our two MONDian implementations. There we 
can see some objects that are still approaching each other in 
OCBMond to have already merged in OCBMond2 indicative 
of an advanced evolutionary stage. 

We quantify the evolutionary stage in Figure[8]where we 
present the matter power spectra at redshift z — 0, z = 2, 
and z — 5. We basically recover the same features as re- 
ported by KG04 (even though they considered a rather dif- 
ferent MOND model with 7(a) = go), namely the lack of a 



Note that A depends both on redshift and cosmology. 
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Figure 7. Density distribution at redshift 2 = (bottom row), 2 = 2 (middle row), and 2 = 5 (upper row) for the 128^ particle ACDM 
(loft column), OCBMond (middle column), and OCBMond2 (right column) simulations. 



distinctive "break" due to the transfer of power from large 
to small scales in the MONDian runs and the (marginally) 
larger amplitude of power on scales k 5C l/iMpc~^ close to 
the fundamental mode. However, we also observe that OCB- 
Mond2 evolved marginally faster on large scale and at later 
times than OCBMond. 

A rather natural question arises about the sites where 
MOND is actually aflective. To shed some hght on this is- 
sue we show in Figure[9]the modulus of the MONDian force 
Sm2 ill i'hs OCBMond2 model at the particle positions at 
redshift 2 = normalized by a^go (which corresponds to the 
argument of the MONDian interpolation function used with 
Eq. Q). We find that particles in low-density regions and 



the vicinity of objects, respectively, are in the MOND regime 
whereas particles residing at the centres of objects (i.e. high- 
density regions) are mostly unaffected by MOND. This may 
result in different (hierarchical?) structure formation scenar- 
ios as matter infall onto objects will most certainly be af- 
fected by the "MONDian environment" of objects. We will 
investigate this in the following Section. 

4.3 Hierarchical Structure Formation 

As already noted in Figure [7] and |8] there appears to be 
a marginally faster evolution of (large-scale) structures in 
OCBMond2. This is further confirmed in Figure [10] where 
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Figure 8. Dark matter power spectrum of all simulations at red- 
shifts z = 5, z = 2, and z = 0. 




Figure 9. Argument of the MONDian interpolation function fJ,{x) 
at redshift z = in the OCBMond2 model, i.e. the modulus of 
the MONDian force | V"3?m| normalized by o?go. The values have 
been interpolated to the particle's positions and are colour-coded 
according to the scale shown on the right. For values above 1 
the particles are treated according to Newtonian physics whereas 
values below 1 indicate the MONDian regime. 



we show the abundance evolution of objects with mass 
M > 10"/i"^ Mq. At around redshift z = 1.5 OCBMond2 
starts to develop more objects. Anyways, the late onset 
and faster evolution of structure f ormation under MQND 
ingeneral was already reported by ISandersI ()200ll ) , iNusserl 
(|2002l ) and KG04. However, here we not only used a dif- 
fernet MOND model but also revised the way of presenting 
the comparison between the two MONDian and the Newto- 
nian runs. As OCBMond and OCBMond2 are simulations of 
the gravitational interactions of baryons whereas the ACDM 
model contains both baryonic and dark matter a fair com- 
parision should correct for that. Or in other words, an ob- 



0.0100 rr 




0.0001 I I I I I 

12 3 4 5 

redshift z 

Figure 10. Redshift evolution of the abundance of objects with 
(baryonic) mass Mi, > 10^^ Mq. Note that the baryonic mass 
for the ACDM model has been estimated by multiplying the dark 
matter mass by the baryonic fraction. 



ject of mass M in one of the MOND runs merely contains 
baryons whereas the corresponding object in ACDM con- 
tains both baryons and dark matter. As the idea of MOND 
is to replace dark matter by a modification to Newtonian dy- 
namics we should "remove" the dark halo from the ACDM 
object when performing a cross-comparison. We therefore 
multiply the ACDM masses by the cosmic baryon fraction 
fb = fib/no and refer to all masses as "baryonic mass" Mh 
even though the interactions of baryons other than gravity 
are not modelled in the simulations presented in this study. 
We note that this is not necessarily correct as galaxies prob- 
ably contain fewer baryons than given by the cosmic baryon 
fraction (e.g., iGottlober fc Yepe3l2007l : lOkamoto et al.l 
,200^ . Our "baryonic" ACDM masses are therefore consid- 
ered an upper limit and are likely to be smaller. 

By revising the mass of ACDM dark matter halos we 
note that the abundance evolution is similar yet difi'erent 
when comparing ACDM against the two MONDian mod- 
els. Similar in a sense that the in KG04 advocated dramatic 
difference at high redshift is not as pronounced anymore. 
However, the actual shape of the Newtonian and MONDian 
curves are though rather different with a more gentle in- 
crease in the number of objects in ACDM (at least at red- 
shifts smaller than z ~ 2.5). 

To gain a better understanding of structure formation 
under MOND we constructed a merger tree for each indi- 
vidual object identified at redshift 2 = 0. We then fitted 
t he universal m ass accretion history formula to the data 
(jWechsler et al.ii2od3 ) 

M{z) = M(,^o)e-"" . (17) 

Applying the half-mass criterion (e.g.. lTormen|[l997l l to the 
fits for obtaining the formation redshift (i.e. the formation 
redshift is implicitly defined via 0.5 = e~"^^j3 we are able to 

* Even though we stored 26 snapshots between redshifts z = 
and z = 5 for each model we prefer to use the fitted mass accretion 
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Figure 11. The "baryonic" masses Mi, q of objects at today's 
rcdshift versus their formation redshift Zf. The soUd Unes are 
best-fit power laws to the scatter plot with the exponent of the 
power law given in the upper right corner of each panel. 




Figure 12. Two-point correlation function at redshift z = 0.0 for 
the 1000 most massive objects. Error bars denote Poissonian error 
bars derived from the number of pairs in the respective distance 
bin. 



OCBMorri 

0CBMord2 




quantify hierarchical structure formation. We expect small 
objects to form first and then successively mer ge to form 
larger and larger entities (e.g.. [Davis et al]|l985l l. This can 
be verfied in Figure [TT] where we plot today's mass against 
the formation redshift of that object. We note a clear hier- 
archical tendency even in both MOND models. To quantify 
differences we fitted a simple power law to the data 



(18) 



where the best fit a is given in the upper right corner of 
the respective panel. We find that especially in OCBMond2 
high-mass objects appear to form at earlier times than in 
ACDM. Note that this does not contradict the dearth of 
massive objects at high redshifts as observed in Figure [TOl 



history as it will provide us with a more precise estimate for the 
formation time. 



Figure 13. The two-point correlation function for objects of dif- 
ferent mass in OCBMond and OCBMond2. The considered mass 
range in the calculation of ^g^i is given in the upper right corner 
of each panel. The error bars are the Poissonian errors again. 



Here we are tracing back individual objects whereas in Fig- 
ure [To] we considered all objects above a given mass at a 
certain redshift; Figure [TT] simply indicates that (some of) 
the high mass objects found at redshift z — must have 
already been in place at a higher redshift than in, for in- 
stance, ACDM. Surprisingly there is also a difference be- 
tween OCBMond and OCBmond2 with the former being 
closer to ACDM. This highlights the relevance of the curl- 
field V X h introduced in Eq. ([SJl for structure formation. 

One of the findings of the study by KG04 was that 
the (formation) sites of MONDian objects are marginally 
more correlated. Can we vindicate that the two-point corre- 
lation function ^gai of MONDian objects (i.e. galaxies and 
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hence the subscript gal) has a larger amplitude (at least) on 
small scales. This can be verified in Figure [12] where we plot 
(.saiix) for the 1000 most massive objects in all three mod- 
els. We confirm the finding of KG04 that OCBMond has 
a (marginally) larger amplitude than ACDM, especially on 
small scales. We further acknowledge that OCBMond2 has 
a substantially larger amplitude! In order to better gauge 
the orgin of the difference between OCBMond and OCB- 
Mond2 we split our objects up into different mass bins (i.e. 
Mb^n = [10^ - 10'°]/i-' M0, Mb^i„ = [10" - 10"]/i"' Mq, 
M^in = [10" - 10^^]/i"^ Mq) and calculate ^gai for each 
mass bin independently. The results can be viewed in Fig- 
ure [TS] We observe that the differences are due to the biased 
formation sites of high-mass objects. We like to mention that 
the number of objects in each mass bin is practically identi- 
cal amongst OCBMond and OCBMond2; only in the range 
10" - lO^^ft-^ Mq there are about 24% more objects in 
OCBMond2 as indicated by Figure 1101 We hence ascribe 
the bias between OCBMond and OCBMond2 as found in 
Figure [12] to the large mass objects and hence the stronger 
evolution of large scales as already seen in Figure [7] and [H] 

4.4 Cross-Correlation of Properties 

Something not considered in the study of KG04 is the cross- 
correlation of individual objects. To this extent we use the 
tool MergerTree that comes with the distribution of the N- 
body code AMIGA and serves the purpose of identifying cor- 
responding objects either in the same simulation at different 
redshifts (and hence the name "MergerTree") or in simula- 
tions of different cosmological models but run with the same 
initial phases for the initial conditions. The cross-correlation 
is done by linking objects that share the most common par- 
ticles. 

We start with the most simple yet still interesting quan- 
tify for our cross-comparison, namely the mass M. In Fig- 
ure Owe show how Mb correlates across ACDM and OCB- 
Mond (upper panel) as well as across OCBMond and OCB- 
Mond2 (lower panel). Remember that we defined Mf, as the 
"baryonic" mass of our halos that coincides with the actual 
total mass in the two MOND models but corresponds to 
A/tot X /i, for ACDM; by the usage of the cosmic baryon frac- 
tion /(, = fii, /Slo in this formula we will have an upper limit 
for the "baryonic" A CDM mass (cf. iCottlober fc YepesI 
l2007l : lOkamoto et al.ll2008i ). Despite the "baryonic" correc- 
tion we still observe a difference between the masses of "iden- 
tical" objects in ACDM and the MONDian runs. This can 
in parts be ascribed to the use of the cosmic baryon frac- 
tion and in parts to the definition of the edge of objects: 
remember that we use A = 340 for ACDM as opposed 
to A = 2200 for MOND. The latter value translates into 
a smaller radius and hence less mass. When applying the 
same density threshold of A = 340 to the MONDian runs 
(not presented here though) the MONDian masses increase 
by approximately 30% bringing them into better (yet not 
perfect) agreement with the ACDM masses. 

The other interesting observation in FigureOis the fact 
that there is a systematical difference between the masses 
of cross-identified objects in OCBMond and OCBMond2. 
While the most massive objects appear to have identical 
masses, there is a clear trend for low-mass OCBMond2 ob- 
jects to be more massive than their OCBMond counterparts. 



ACDM vs. OCBMond 




Figure 14. Ratio of masses for cross-identified objects at redshift 
2 = 0. The histograms represent the mean ratio in the respective 
bin. Note that the ACDM masses have been lowered by the bary- 
onic fraction to be comparable to the MONDian values. 



This difference fiattens to values closer to unity when con- 
sidering objects at higher redshift. This phenomenon can 
therefore be ascribed and explained by the even stronger 
evolution of the OCBMond2 simulation as already noticed 
in Figure [7] and Figure 1101 respectively. As pointed out by 
Figure [9] MOND is most affective in the outer regions of 
objects and hence leaving its imprint via differing infall pat- 
terns of material. We therefore ascribed the differences to 
variation in the mass accretion histories as confirmed by 
Figure 1111 

We close this Section with a brief investigation of the 
shape of o bjects as defined b y, for instance, the triaxiality 
parameter (|Franx et al.lll99ll ) 



where a > b > c are the principal axes of the moment of 
intertia tensor we note that there are hardly any differences 
at all. This can be verified in Figure [15] where we plot for 
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Figure 15. Ratio of triaxialities for cross-identified objects at 
redshift z = 0. The histograms represent the mean ratio in the 
respective bin. 



the same set of cross-identified halos as already used in Fig- 
ure [Utile ratio of the respective triaxialities as a function 
of mass again. This figure highlights that despite variations 
in the mass growth of objects at least their triaxialities are 
unaffected. 



4.5 The curl field 

As we have just seen there are subtle yet noticable differ- 
ences when using our novel Poisson solver for the MONDian 
equation Q as oppossed to the simplified prescription of 
inverting Eq. under the assumption of C = 0. In that 
regards it appears mandatory to study the origin of these 
devations. We therefore present in this Section a detailed 
investigation of the curl-field V x h and its effects upon 
structure formation. 

Given the solution of the Eq. ([2]) and Eq. Q the curl- 
field is readily calculated 




^lgM2 




gN = ji9M 



Figure 16. Sketch illustrating the vectors of Eq. I I20I I and the 

angles between them, gjy is the Newtonian force vector, gj^j the 
one used with the OCBMond model, and gj(/2 the solution of our 
new solver responsbile to the evolution of the OCBMond2 model. 



V X h = 



(20) 



Note that setting V x h = corresponds to the OCBMond 
model. 

In order to clarify all definitions and terminology to be 
used from now on we show in Figure [16] a sketch of all the 
vectors gj^ = -V$jv, gj,^ = -V$m, and g^.^j = -V$m2 
in play. The assumptions used with the OCBMond model 
imply that the forces there and the Newtonian forces, gj^j 
and gjY, are parallel whereas they both can have an angle a 
with the OCBMond2 force vector gj^.i2 ■ This angle is due to 
the influence of the curl-field and can be described by the 
angle (3 (cf. Eq. CT). 

Interpreting the curl-field as a modification to gjv on 
the right-hand side of Eq. ([SJ the question arises about the 
actual adjustment to (the implicitly defined) gjjjj stemming 
from this additional "source" term. We therefore normalize 
C by the Newtonian acceleration gjy and study its spatial 
dependence. This can be viewed in Figure [17] where we show 
|C|/|gjY| (lower panel) in a slice of thickness 125/i~^ kpc in 
the OCBMond2 simulation as calculated on a 256^ grid. We 
do observe a marginal correlation of the "correction to the 
source" with the actual matter density shown in the upper 
panel. We further note that the modulus of the curl-field 
is actually always smaller than the Newtonian force, how- 
ever, its spatial variation is though highly non-trivial. In 
some (high-density) regions there appear to be some kind of 
"shocks" where the curl-field changes dramatically changes 
its value. The only conclusion we can draw from this qualita- 
tive analysis is that the curl-field has an infiuence primarily 
in low-density regions: its value is close to the Newtonian 
force in the void regions. Finally we like to remark on the 
"diamond shapes" seen in Figure [17] they are an artifact of 
the periodic boundary conditions and stem from the peri- 
odic images of the largest object /void. These attracting im- 
ages create the observed symmetry leading to the diamond 
shapes. 

In order to make a more quantitative analysis of the 
curl-field, we use the snapshots of the simulation OCB- 
Mond2 and calculate the Newtonian as well as the two 
MONDian forces at various redshifts. These forces are then 
extrapolated to the particle positions in the same way as 
done during the time integration of Eq. ([TJ. We are aware 
that this extrapolation will bias our results towards high- 
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Figure 17. Slice of tliickness 125/i-l kpc of tlie 0CBMond2 sim- 
ulation at redshift z = 0. The slice is chosen to contain a strong 
density peak (i.e. it is centred about the position of the most 
massive object found in th simulation). The upper panel shows 
the density as obtained on a regular 256^ grid. The same grid 
is used to calculate the modulus of the curl field normalized by 
the modulus of the Newtonian force, i.e. |V X h|/|gjY|, which is 
presented in the lower panel. 



density regions (where the particles actually reside) and 
ignore the efFects-of-interest in low-density regions, respec- 
tively. However, our integration scheme is based upon the 
idea of sampling phase-space with particles and hence the 
forces at the particle positions are relevant for the time inte- 
gration of Eq. lU and hence the evolution of the simulation. 

In Figure [18] we now present a quantitative comparison 
of the curl field, curl = V x h, and the Newtonian forces 
Siv by plotting the probability distribution of the angle be- 
tween the two vectors (upper panel) as well as the relative 
difference of the modulus (lower panel). First, we note that 
there is practically no evolution with redshift. Second, the 
curl-field is well aligned with the Newtonian force though it 
may also point in the opposite direction! In order to empha- 
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Figure 18. Probability distribution of the angle between the curl- 
field V X h and the Newtonian force gjy (upper panel) and the 
fractional difference between these two values (lower panel) at 
various redshifts. 



size the skewness of the distribution about cos (curl, g^) = 
we show the modulus of it. The upper curves thereby corre- 
spond to cos (curl, g^) < 0, i.e. it is more likely that the curl 
points in the opposite direction. However, the lower panel of 
Figure [18] indicates that the actual change to the "source" 
term for the implicitly defined gj^^j (cf- Eq. ((5]l) induced by 
the curl is rather small. Here we show the fractional differ- 
ence between the moduli of the curl and the Newtonian force 
vector. This distribution peaks at approximately —0.9 and 
hence a ~10% correction to g^ in Eq. ([6]l; this modification 
is again independent of redshift. So, the net effect of the 
curl field is to reduce the "source" on the right hand side of 
Eq. ((6]l . We expect this to translate into a reduction of gj^^j 



with respect to gJ^^, too. 
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Figure 19. Probability distribution of tlie angle between the 
MONDian forces gj^j and g]^,J2 (upper panel) and the fractional 
difference between these two values (lower panel) at various red- 
shifts. 
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Figure 20. Projection of a sub-box of side length 1.5h~ Mpc 
at rcdshift 2 = 0. The left plot shows a colour-representation 
of cos (gj\/,gjvf2) whereas the right plot shows the density field. 
Both values are evaluated at the particles' positions. Note that 
both figures are centred on the same point. 



towards values of stronger gj^j and lower gj^j2' respectively. 
This skewness is a manifestation of the result obsvered in 
Figure [181 namely that the curl preferentially decreases the 
"source" in Eq. ©. 

We further note in Figure [19] that there is a marginal 
redshift evolution of the peak in the distribution of the frac- 
tional differences of the two MONDian forces. The peak 
gradually migrates away from at redshift z = 5 towards 
more negative values of approximately —0.09 at z = 0. While 
it appears counter-intuative to explain the advanced struc- 
ture formation in OCBMond2 found in previous Sections 
(cf. Figure [7] Figure [S] Figure 1101 etc.) with the finding 
that the forces are smaller in that model, we rather believe 
that the shift in the peak can be held responsible for it. 
This shift towards negative values translates into stronger 
OCBMondS forces and hence structure formation at an ac- 
celerated speed. 

As an illustrative example of the misalignment between 
the two MONDian forces (and its relation to the underlying 
density field) we show in Figure [20] both the (colour-coded) 
cosine of the angle between gj^j and g^^j (l^ft panel) along- 
side the density at the particles' positions at redshift z = 0. 
We note that there is no apparant correlation of the mis- 
alignment with density. 



But because of the implicit definition and the vector- 
nature of the quantities involved in the calculation of gj^f2 it 
is difficult to make predictions for the change in gJ^.J2 induced 
by the modifications to the "source" in Eq. ((Gjl . We therefore 
plot in Figure[T9]the analogous quantities as in FiEure [T8l this 
time for the comparison of gj^j and gj^j2' ^^is should reveal 
the effects of the curl-field directly. As the angle distribution 
is no longer symmetric we expand it over the whole range 
from cos (gj^j^, g^jj) ^ The upper panel (showing 

the angle distribution) clearly indicates that both MONDian 
forces are well aligned (note the logarithmic scale on the 
2/-axis). The lower panel (showing the fractional difference 
between gj^^ and gM2) proofs what we already expected: 
the distribution is centered about but shows a skewness 



5 SUMMARY & CONCLUSIONS 

We presented a novel solver for the analogue to Poisson's 
equation taking into account the effects of modified New- 
tonian dynamics (MOND). This equation is a highly non- 
linear partial differential equation for which standard solvers 
based upon Fourier transformation techniques and/or tree 
structures are not applicable anymore; on e has to de- 
fer to multi-grid relaxation techniq ues (e.g., iBrandtl Il977l : 
[Press et al. I I 19921 : iKnebe et al.|[200ll ) in order to numerically 
solve it. 

The major part of this paper hence deals with the nec- 
essary (non-trivial) adaptions to the existing multi-grid re- 
laxation solver AMI GA (successor to MLAPM introduced by 
iKnebe et al.l (|200ll )). We show that the accuracy of our 
MOND solver is at a credible level for a static problem with 
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known MONDian solution and that it reproduces correct 
results in the Newtonian limit for cosmological simulations. 

In the case that we don't want to adhere to any (und- 
justifyable and hence ad-hoc) assumptions as, for instance, 
done by KG04, we are facing the problem that there is still 
a lot of freedom from the covariant point of view to describe 
the original phenomenological MOND theory. To be able 
to actually perform MONDian cosmological simulations we 
decided to choose one particular model for MOND and to 
leave for later studies the analysis of the effects that dif- 
ferents theories producces on the structure formation. One 
of the suppositions KG04 made relates to the relation be- 
tween the MONDian and the Newtonian force vector. They 
use the rather simple relation as given by setting C = in 
Eq. (|6]) making the right-hand side simply the Newtonian 
force. This implicit definition for ^j^j could be inverted and 
used instead of gjv when updating the particles' velocities 
during the time integration of Eq. ([T}. As noted in Eq. ([6]) 
they neglected the so-called MONDian curl-field C = V x h. 
Others have shown that this field de creases like 0(r~^) and 
vanishes for any kind of symmetry (|Bekenstein fc MilgromI 
1 19841 '). it yet remains unclear whether it will leave any im- 
print on inhomogeneous strucure formation as found in cos- 
mological simulations. 

In the second part of this paper we therefore presented 
a series of cosmological simulations set out to quantify both 
structure formation under MOND as well as the effects of 
the curl-field. Surprisingly, we found that our results are con- 
sistent with KG04 even if we take in account that they use a 
different MOND model and initial conditions generated with 
standard linear theory in a MONDian regimen. We revised 
some of their findings (mainly the discrepancy between the 
abundance of galaxies at galaxies at high redshift between 
ACDM and OCBMond) and noted that the curl-field leaves 
marginal yet noticable effects. The major result of this study 
though is that the curl-field appears to drive structure for- 
mation! Our results obtained with the new solver (and hence 
including the effects of the curl-field) can be summarized as 
follows 

• the curl-field drives structure formation, 

• the curl-field leads to more objects at 2 = where 

• cross-identified objects are more massive in OCBMond2 
than in OCBMond, and 

• the OCBMond2 model shows a stronger two-point cor- 
relation function. 

We acknowledge that there are still a lot of puzzling 
results to be investigated in greater detail. However, we 
postpone this to future studies. The aim of this paper was 
primarily to describe the novel gravity solver that is freely 
available for download 



^ AMIGA can bo downloaded from the following web page 
http://Mww.aip.de/People/aknebe/AMIGA. The new MOND 
solver can be switched on via the compilation flag -DMQND2. The 
code is able to run cosmological simulations, ie with expansion 
and periodic boundary conditions and to work also in isolation, 
wich implies no expansion and fixed boundary conditions. It can 
be also use without temporal evolution in order to get 3D po- 
tentials from density distributions given by particles or analitic 
formulas. Please contact the authors for more details. 
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